#ifndef __QUADRATURE_H__
#define __QUADRATURE_H__

#include <cstdlib>
#include <iostream>
#include <stdexcept>
#include <valarray>
#include <string>
#include "Point.h"

/* ------------------------- class Quadrature begin --------------------------- */
template <size_t DIM>
class Quadrature
{
protected:
	size_t accuracy;
	std::valarray<Point<DIM, double> > quad_points;
    std::valarray<double> quad_weights;

public:
	Quadrature() = default;
	virtual ~Quadrature() = default;
	size_t GetAccuracy() const;
	size_t GetNQuadPoints() const;
	const std::valarray<Point<DIM, double> >& GetQuadPoints() const;
	const std::valarray<double>& GetQuadWeights() const;
	virtual void ReadQuadInfo(const std::string& file_name, size_t acc) = 0;
};

template<size_t DIM>
size_t Quadrature<DIM>::GetAccuracy() const
{
	return accuracy;
}

template<size_t DIM>
size_t Quadrature<DIM>::GetNQuadPoints() const
{
	size_t n_q_pnts = quad_points.size();
    size_t n_q_wei = quad_weights.size();
	if (n_q_pnts != n_q_wei) {
		throw std::length_error("number of quad_points != number of quad_weights!");
		std::exit(-1);
	}
	return n_q_pnts;
}

template <size_t DIM>
const std::valarray<Point<DIM, double> >& Quadrature<DIM>::GetQuadPoints() const
{
	return quad_points;
}

template<size_t DIM>
const std::valarray<double>& Quadrature<DIM>::GetQuadWeights() const
{
	return quad_weights;
}
/* ------------------------- class Quadrature end --------------------------- */



/* ------------------------- class TriangularQuad begin --------------------------- */
class TriangularQuad : public Quadrature<2>
{
public:
	TriangularQuad() = default;
	virtual ~TriangularQuad() = default;
	virtual void ReadQuadInfo(const std::string& file_name, size_t acc);
};
/* ------------------------- class TriangularQuad end --------------------------- */


#endif
